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Abstract. The purpose of this paper is to improve the modeling of the mixing of chemical elements that occurs in stellar 
radiation zones. In addition to the classical rotational mixing considered in our previous paper, which results of the combined 

■ action of the thermally-driven meridional circulation and of the turbulence generated by the shear of differential rotation, we 

include here the effect of an axisymmetric magnetic field in a self-consistent way. We treat the advection of the field by the 
7-H ■ meridional circulation, its Ohmic diffusion, and the production of its toroidal component through the shear of differential 

' rotation. The Lorentz force is assumed not to exceed the centrifugal force; it acts on the baroclinic balance and therefore on the 

meridional flow, and it has a strong impact on the transport of angular momentum. All variables and governing equations are 
expanded in spherical or spherical vectorial functions, to arbitrary order: this yields a system of partial differential equations in 
time and in the radial coordinate, which is ready to be implemented in a stellar structure code. 
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The need for improved stellar structure models, which go beyond the so-called standard model, is now well recognized. The 



1. Motivation 

Of 
i 

2 ' radiation zones can no longer be treated as motionless regions, where no mixing occurs, but one has to include physical processes 
\ which lead to the transport of matter and angular momentum. So far the focus has been on the effect of rotation, which generates 
d • a thermally driven meridional circulation; by advecting angular momentum, that circulation renders the rotation non-uniform and 
| prone to shear turbulence. This rotational mixing has been formulated by Zahn (1992) and Maeder & Zahn (1998), assuming that 
• *~j ■ the angular velocity varies much less in the horizontal than in the vertical direction (shellular rotation). For rapidly rotating stars, 
r> ' models built along these lines are in much better agreement than standard models; this was first shown by Talon et al. (1997), 
^ ■ and confirmed by Maeder and Meynet (2000), Maeder and Meynet (2001), who in addition refined the description for the mass 
loss by taking into account the latitudinal variation of the wind. 

However, these models fail to correctly reproduce the flat rotation profile observed in the solar radiative interior (Matias & 
Zahn 1997). In solar-type stars, at least once they have been spun down through their magnetized wind, other physical processes 
are therefore responsible for the transport of angular momentum, and the most plausible candidates are now being investigated. 
One is the transport by internal gravity waves emitted at the base of the convection zones, associated with turbulence which 
smoothes the differential rotation; this mechanism was described by Kumar et al. (1999), Talon et al. (2002) and applied by Talon 
and Charbonnel (2003, 2004) to stars of various types. The other possibility is magnetic torquing, as advocated already by Mestel 
(1953): he showed that even a very modest field could enforce rigid rotation. In the present paper we shall address the latter 
effect, and examine how rotational mixing is affected by an axisymmetric field. 

Moss (1974) was the first to map the thermally-driven meridional circulation in a star containing an axisymmetric poloidal 
field; when rotation was added, it was assumed constant and uniform. Various generalizations of this work followed, by Mestel 
and Moss (1977) and by Moss (1977, 1984, 1985, 1987). The a priori assumption of rigid rotation was relaxed by Tassoul 
and Tassoul (1986), who claimed that no magnetic field could render the rotation nearly uniform in presence of meridional 
circulation; however their treatment failed to correctly represent the generation of toroidal field through the shearing of poloidal 
field, which provides the major feed-back in this problem. This was pointed out by Mestel et al. (1988), who carried out numerical 
simulations that illustrated how a weak magnetic field could enforce nearly uniform rotation within a few Alfven times; for sake 
of simplification, both the poloidal field and the meridional circulation were taken as given, and constant in time. 
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Recently Garaud (2002) treated the problem in a fully consistent way, and she applied it to the description of the solar 
tachocline: she let the field adjust through Ohmic diffusion and be advected by the meridional circulation. However, due to 
numerical limitations, her circulation was driven mainly through Ekman pumping, rather than through thermal diffusion; also, 
only stationary solutions were constructed. 

The purpose of the present paper is to model the time-dependent rotational mixing in presence of an axisymmetric magnetic 
field, thus extending the results obtained by Mathis & Zahn (2004; from now on referred to as Paper I) in the purely hydrody- 
namical case. To achieve this goal: 

- we allow for a general axisymmetric magnetic field, including all its effects, 

- we expand all variables in spherical harmonics, formally to unlimited order, 

- we treat explicitely the departures from shellular rotation, in the linear approximation, in order to capture the tachocline(s), 

- we filter out the short times, and keep only the time derivative in the heat equation and in that describing the transport of 
angular momentum. 1 

In a forthcoming paper, we shall generalize the present work to a non-axisymmetric magnetic field. 

2. Assumptions and expansions 

We focus here our attention to magnetic fields of moderate strength, whose Alfven speed Va does not exceed the rotation velocity 
RQ, but is larger than the meridional circulation velocity Vm- Thus the field will not affect much the hydrostatic balance in the 
meridional plane, but it may play a major role in the transport of angular momentum, by tending to render the angular velocity 
O constant along the field lines of the poloidal field (Ferraro 1937). We also assume that the field does not vary on a timescale 
which is shorter than the Alfven time R/Va', thus we shall deal with fields of primordial origin, and exclude possible dynamo 
fields produced in the adjacent convection zones, whose penetration would be damped anyway by skin effect (Garaud 1999). 

Our treatment is not able to address the question whether the field is stable or not, because it is axisymmetric, and because we 
do not resolve the Alfven waves. We will have to rely on our intuition, guided by 3D simulations such as recently performed by 
Braithwaite and Spruit (2004), to choose the magnetic configurations we implement in our code. As for instabilities which may 
lead to a turbulent state, for instance the magneto-rotational instability when the equator rotates slower than the poles (Balbus & 
Hawley 1994), their transport properties will have to be accounted for by a suitable parametrization, as was done previously in 
the hydrodynamical case. 

As in Paper I, we assume that the rotating star is only weakly two-dimensional, and this for two reasons. The first is that 
the rotation rate and the magnetic field, which is taken axisymmetric, are sufficiently moderate to allow the centrifugal and the 
Lorentz forces to be considered as perturbations compared to gravity. The second reason rests on the less justified hypothesis that 
the shear instabilities due to the differential rotation give rise to turbulent motions which are strongly anisotropic due to the stable 
stratification, with much stronger transport in the horizontal directions than in the vertical. In the radiation zones, we expect it to 
smooth out the horizontal variations of angular velocity and of chemical composition, a property we shall use to discard certain 
non-linear terms. Let us emphasize also that the influence of the magnetic field on such turbulence is not taken into account, as 
seems to be allowed by the condition Va < RQ- 

Hence, we consider an axisymmetric star, and assume that the horizontal variations of all quantities are small and smooth 
enough to allow their linearization and their expansion in a modest number of spherical harmonics. As reference surface, we 
choose either the sphere or the isobar, and write all scalar quantities either as: 




(1) 



or 




(2) 



where Xi(r) and Xi(r) are related by (cf. Paper I, §2) 




(3) 



with r being here the mean radius of an isobar. 



1 All waves are thus filtered out (acoustic, gravity, inertial, Alfven); the equation for the transport of angular momentum (eg. 1441 below) 
still formally allows for mixed Alfven-inertial waves, but the timestep of the evolutionary calculation (which uses an implicit scheme) will in 
general exceed the typical wave travel time, and suppress these waves. 
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Likewise, we expand all axisymmetric vector fields using the method outlined by Rieutord (1987): 

u(r, 9) = £ ("oW*z W) + v*,(r)S°(0) + w l (r)T^(9)} . (4) 

1=0 

The orthogonal axisymmetric vectorial spherical harmonics /?? (#)> S9 (9), J? (0) are defined by: 

R ( }(0) = Y°(6)? r , S°(0) = V S Y°(9) and T°(6) — A (5) 
where Y? (9) is the classical spherical harmonic with m = (cf. Edmonds 1974) and is the horizontal gradient 

Vs = 7 g d e +^-^—8^. (6) 
sin# 

Their detailed properties are given in Appendix A. 

Using IA.H and ( IA.6> . « may be written explicitly on the usual vector basis in spherical coordinates: 



" = Z ^ K (r) P/ (cos + v '> (r) 5flP/ (cos 9)Fe ~ w 'v (r) 5flP/ (cos 

1=0 

CO 

= A$k° (r) % + J] < [4 (r) P, (cos 0) e" r - v(, (r) P> (cos 9)T e + w l (r) P> (cos 0)T V ] ; (7) 
z=i 

~e r , Tg and<T^ are the unit- vectors repectively in the r, 9 and <p directions. 

This expansion of the vector fields allows us to separate explicitly the spatial coordinates (r, 9) in the vectorial partial differ- 
ential equations which govern the problem: using this decomposition, we cast the problem into partial differential equations in t 
and r only. This point is very important numerically, because the existing stellar evolution codes in which we have to introduce 
the transport equations are all 1 -dimensional, and also because we need to achieve much higher accuracy in the radial direction 
than in the horizontal. 



3. Axisymmetric magnetic field 

We start with expansions and equations involving the magnetic field. 



3.1. Definitions 

The magnetic field, being divergenceless, is conveniently split in its poloidal and toroidal parts: 
B (r, 9) = B P (r, 9) + B T (r, 9) = V A V A (fr (r, 0)T r ) + V A (fr (r, 9)T r ) , 



(8) 



where %p and £7- are respectively the poloidal and the toroidal magnetic stream-functions. Using the classical method introduced 
by Bullard and Gellman (1954) for the spectral treatment of the geodynamo problem (see also James 1973, 1974 and Serebrianaya 
1988), we get the following expansion in radial functions and in spherical harmonics which are defined in §A.1.1.: 



(r, 6) = J] ^ (r) Y° (9) and & (r, 0) = J^o W *? ( e ) ■ 



(9) 



/=I 



Since R.9 (9) = if (6) e,-, as we have seen in the previous section where we have defined R.9 (9) , S' 1 (9) , T® {9), the field may be 
written as 

CO OO 

B, = 2 VAVA (^?®) and B T =2 VA (^ (0) )- (10) 

/=i z=i 

Next, we develop the curl operator as explained in the appendix (eqs. jA.32> and ( IA.33t ) to obtain the following expansions for 
the magnetic field: 



B (r, 9) = Bp (r, 9) + B T (r, 9)=^- 



mvfe 



Wo 



1=] 



X'o 



(11) 



We proceed likewise for the current density /, which in the framework of magnetohydrodynamics is related with B through 
the Maxwell- Ampere equation, neglecting the displacement current: 



j(r,0) = — VAfi(r,0, 



(12) 
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/io being the magnetic permeability of vacuum. 

Therefore, using expression for the field, we get: 
j - —V A B - — [VaVaVa (^pe r ) + V A V A (£re"r)] , 

MO MO 



(13) 



which, using the relations JA.33> and ( IA.36t . we can project again on R® (0), S? (ff) and Ty (0), and split into its poloidal and 
toroidal parts: 



1(1 + 1)4 



1 1 j 

j P (r, ff) — — V A B T (r, ff) — — V A V A = — Y 

Mo MO Mo 

1 1 1 °° f 

/V (r, 0) = —V A B/> (r, 8) = — V A V A V A (£ P F r ) = V i 

Mo Mo Mo ^ 



( c l Y 



- d rXo 

r 



S a ,(0) 



r?(e) 



with A/ being the laplacian operator 
A, = d r>r + -d r — . 



(14) 



(15) 



(16) 



We are now ready to examine the aspects of transport and mixing in the radiation zones of rotating stars that are related to 
the presence of an axisymmetric magnetic field. To achieve this, we shall proceed in three steps 

- first, we introduce the Lorentz force; 

- next, we derive the transport equation for B, which is the classical induction equation; 

- finally, we take into account the energy losses due to Ohmic heating. 



3.2. Lorentz force 

The Lorentz force fx plays a crucial role in a rotating star, since it acts to render the angular velocity constant on the field lines 
of the poloidal field Bp. In this section, we establish the formal expression of 9x i n terms of vectorial spherical harmonics. 
Starting from the definition of the Lorentz force: 



T £ (r,6)=jAB 



— (V A B) 

y"0 



AB, 



(17) 



we replace B by its expansions ( II It . and use the algebra ruling the vector product of two general axisymmetric vectors in 
Appendix A, whereby we obtain the following formal projection: 

CO 

Tz (r, 8) = Xf x ;o(r)Tr + J] [Xf x j(r)Pi(cos 8Je r + Jf fj .^r)P}(coa ffjei, + Zt z j(r)P)(cos 8)^) 



i=\ 



Xr £ -,o(r) 



*o«9 + £ 
1=1 



Xr z j(r) 



Nf 



R°(0) + 



Nf 



■w + 



Zr £ ;l( r ) 



Nf 



(18) 



In Appendix B, we give the explicit expressions for the radial functions X^ £ -i (r), ifj^ j ( r X Zr £ j ( r ) i n terms of the magnetic 
stream-functions £ Q and^,, the normalization coefficient Nf being given in (IA.2> . 

As we have done previously for the magnetic field B and for the current /, we split the magnetic force into its poloidal and 
toroidal parts: 

OO 

T L ,P = j T A Bp + j P A Bp — Xr^r + {XTfjPfcoS SWr + ^TzjPj (COS 



1=1 



Xtt-.o 



i=i 



x 



Nf 



R,(0) + 



Nf 



T £J =j P ABp = Y J Zt £ jPi( cos0 K = Yj 



i=] 



Zr, 



Nf 



(19) 
(20) 
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As it was underlined by Mestel et al. (1988) and as we shall see in §4, §5 and §6, this decomposition is very useful: the 
poloidal Lorentz force Tj^ p operates on the hydrostatic balance and thus contributes to the thermal imbalance, while the toroidal 
component Tcj ac ts on the transport of angular momentum through its torque. 

The expression of Ti_ in the case where only the I — { 1 , 2, 3} terms are kept in the expansion of B can be found in Appendix E. 
One should note that if the expansion of B is terminated at mode Z max = N,„ in (II lb . then the expansion of Tz involves terms up 
to / m ax = 2N m , due to selection rules (see Appendix B). 



3.3. Induction equation 

We now turn to the evolution of the magnetic field in stellar radiation zones, which is governed by the induction equation: 

d,B = V A (V AB) - V A (IfolloV AB), (21) 



where we have allowed for an anisotropic eddy diffusivity \\r]\\. We recall (cf. Paper I) that the macroscopic velocity field V is the 
sum of a zonal flow = r sin 6 Q. (r, 0)~e^ and of a meridional flow 11 (r, 0): 



V(r,9) = r sin Q.(r, 0)7^ + U(r,0). 



(22) 



The latter can be split into a spherically symmetric part, which represents the contractions and dilatations of the star during its 
evolution, plus the thermally-driven circulation 



<U{r,e) = r-e r + f U M (r,e), 

the meridional flow being expanded in spherical functions: 



/>o 



_ dP,(cos 0)_ 

U,(r)P, (cos 0) e r + V,(r) '\ ' e g 

d6 



(23) 



(24) 



The radial functions t//(r) and V/(r) are related by the continuity equation, i.e. V ■ (pUlu) — in the anelastic approximation: 
1 d 



1(1 + l)prdr 

We introduce the expanded form of the velocity field in the induction equation d2 It : 
d,B - V A (fT r A B) = V A [fay + <U M ) A B] - V A (| Ml ® V A B) , 



(25) 



(26) 



and as before we project this equation on the vectorial spherical harmonics R® (6), S® (9), 7*} (0). The time derivative on the 
left-hand side of J26l > is readily derived (cf.ll 1>: 



Kl+D d -$ 



-d,A 



2f(fl)L 



(27) 



Then, using the algebra related to the /?? (0), S? (0), 7*} (0) and the identity JA.32> . the following expansion of the contrac- 
tions/dilatations term is obtained: 



V A (rTr AB) = ~Y J l [1(1 + D^34] *° (0) + -d r (rd^ ) S° (0) + -d r (r X {) 



t(0)\. 



(28) 



In the same way, we get for the dissipation term: 



VA(|to||®VAB) = - J] 



1=1 



r 



1(1+1) 



R°(0) + 



-d r 



sU0) 



-d r (r] h d r x' Q )-r] v l(l+l) 



(29) 



where we recall that A/ is the laplacian operator: A; = d rj + (2/r)d r - 1(1 + l)/r 2 . 

The task is more difficult with the non-linear expressions describing the stretching and advection of the magnetic field: 
V A [Ktp A fij + V A [Km A B]. The first term corresponds to the generation of toroidal field trough the shearing of the poloidal 
component by the differential rotation, while the second one represents the advection of the field by the meridional circulation. 
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The first step is to expand {14 v + 14 M ) on the R° {0), S° (0) and T° (0): 

CO 

K v + 14 M = J] [u' Q (r) *° (0) + y(, (r) S° (0) + w (r) 7° (0)} . 



From J24l >. we get immediately for the meridional part: 

l f . U,(r) , V,(r) 
u o W = -Jfi- and v o to = 



Nf 



Next, we expand the rotation law as in Paper I: 

CO 

D. (r, 6») = Q(r) + Q (r, 0) = Q.{r) + Y Q, (r) g/ (0) = ^ Q,* (r) P/ (cos 0) , 



(30) 



(31) 



(32) 



1=0 



where Q (r)= J () Q (r, 0) sin J 0d0/ J () sin 0d0 the horizontal functions are given by Qi (0)-Pi (cos 6*) - //, with 
/, = JT P, (cos 0) sin 3 0d0/£ sin 3 0d0 = 5 W - \d k2 . Thus, we get: 

f Qq (r) = Q (r) - 2;>o to // = flo (r) + ±Q 2 (r) 
\ Q; (r) = Q; (r) for 2 > 0. 

Using the identity ( IA.9> . we obtain the zonal flow 14^ = Yii>o w b ( r ) (0)' where w' Q (r) is given by 



(33) 



(34) 



The expression of the wi when we keep only the two first terms of the expansion (I32> . O and O2 is given in JE.3i . 

It remains to project the vector product of the two axisymmetric vectors, (14^ + 14m) and B, on the vector harmonics R® {0), 
(0) and J? (0)- This task is accomplished in §A.2.3.; applying it to the advection term, we get: 

00 

(tt v + 14 M ) Afi = ^ {^Ad;z(r)^(cos 0)T r + ^Ad;i{r)P}{cos 6)eg + ZAd;i(r)Pj(cos 0Je^} 



-S- 



Nf 



R°{0) + 



Nf 



S°M + 



Nf 



(35) 



where the radial functions Xa&j (r), J/Ad;/ to an d -2 Ad;/ (O, which are explicit functions of Q, Q/, Ui, A and K,, are given in 
Appendix C. Finally, using once again (IA.32> . we obtain: 



V A [(l4 v +14 m )aB\ = Yj 



Nf r 



R l !(0) 



l / Xaaj l 
h -i 

r r 



Jfi\— - + - d r( r ^Ad;l) 



jja-d r (rZw) 
t(0)\- 



(36) 



We are now ready to put in their final form the one-dimensional advection/diffusion equations respectively for the poloidal 
A and for the toroidal magnetic stream-functions x' , by collecting the Rf {0) and the (0) components of J27K d28i . ( I29l l and 
06}: 



"dT = Jp Zm + %rA ' 



dy 1 

+ 5, (r)V = — [X Ad; , + 3,. (rJ/ Ad; ,)] + 
df TV, 



5 r (%3 r y )-77 v /(/+l)^ 



(37) 



(38) 



Note that these equations involve the Lagrangian time-derivative, which makes them suitable for their implementation in stellar 
evolution codes. They are the equivalent in our formalism of the two classical equations for Bp = V A A, A — AT V being the 
potential vector, and Bj - B-fe^. 



d,A + -14 M ■ V(jA) = rj{v 2 A - — j and d,B T + s14 M • v (-y) + B r v ' U P = sB P ■ VQ + j 1 (v 2 B t - ^j, 



(39) 
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where s — r sin and rj is taken isotropic and uniform (Campbell 1997, Mestel 1999). 

Finally, we note that in the ideal case, in the absence of shear, meridional circulation and Ohmic diffusion, fl37i and fl38i 
reduce to: 



— — = and — 
df df 



4" 
r 2 p) 



= and thus to — ( r 2 B r ) = and ~ pUo. (40) 

at V / dt \ rn I 



df V ' df \rp 

These equations express the Lagrangian flux conservation of respectively 

- the magnetic field through the sphere of radius r, as the star expands or contracts, 

- the toroidal field in a homothetic contraction or expansion where the density varies as r~ 3 (cf. Cowling 1957). 

In Appendix E.l. we give the explicit expressions of J37i and J38l > for the dipole, the quadrupole and the octupole (/ = {1, 2, 3}), 
retaining the lowest order terms in the rotation law: Q(r, 6) = Q. (r) + D.2 (r) Q2 (0), and the associated meridional circulation. 

3.4. Ohmic heating 

The last point we shall consider for completeness is Ohmic heating, which contributes somewhat to the thermal bal- 
ance/imbalance. In the case of an anisotropic eddy-magnetic diffusivity, its input rate is given by 

J(r, 0) = — [\\tj\\ ® (V A B)] ■ (V A B) , (41) 
Ho 

which reduces to the classical result J = j 2 /cr, with cr — 1/ (porj) being the conductivity, when the diffusivity is isotropic. 
The scalar product of two axisymmetric vectors is treated in the §A.2.3.; applying the result to J~, we get: 

J (r, 0) = J] J, (r)P, (cos 0), (42) 

/>o 

where the radial functions STi(r), which are functions of the magnetic stream-functions A and K,, are given in Appendix D, and 
in Appendix E for the special case where only the I = {1, 2, 3} modes are kept. Note again here that if the expansion of B is cut 
at mode l m . dx = N m in il It . then the expansion of J' involves terms up to / max = 2N m . 

We are now ready to introduce the Lorentz force in the equations governing the transport of momentum and of heat. 

4. Transport of angular momentum 

We start from the momentum equation 

p [d,V + (V ■ V) V] = -VP-pV0 + V-IMI + Ti, (43) 

where p is the density, the gravitational potential, ||t|| the turbulent stresses and fx me Lorentz force. We insert in it the 
expression for the macroscopic velocity V given above in J22I I23L and take its azimuthal component, which yields an advec- 
tion/diffusion equation for the angular momentum density: 



p^-(r 2 sm 2 0Cl) + V -(pr 2 sin 2 0m( M ) = ^r^-d r (pv v r 4 d r Q) + —dg(pv h sin 3 0d g Q) + r rr (r,8) . (44) 

df V ' V ) r 2. \ I sm Q \ I - 



Note that here again we have introduced the Lagrangian time derivative, making use of the anelastic continuity equation. This 
equation is of course identical to that given in Paper I (eq. 14), except for the magnetic torque Y<f L (r, 0) = r sin 0e<p ■ fx ( r > 
As in Zahn (1992), we assume that the effect of the turbulent stresses on the large scale flow is adequately described by an 
anisotropic eddy viscosity, whose components are v v and Vh respectively in the vertical and horizontal directions. In the absence 
of circulation, turbulence and magnetic field, we retrieve the Lagrangian conservation of angular momentum. 

It remains to project this equation on spherical harmonics; the expansions of the meridional circulation and the angular 
velocity were established in Paper I and they have been recalled above in ( 1241 1321 . We proceed likewise for the magnetic torque: 



Y n (r, 0)=r sin Tx,<p = r sin V Z>F L -\(r)P\ (cos 0) , (45) 



where Tx,ip is the azimuthal component of Tx obtained from ( fT^l . Then, using the following property of the associated Legendre 
functions: 



d d 

— P/ (cos 6») = sin6»— . 
dd dp 



P] (cos 6») = -^Pi (cos 60 = sin 0^-Pi (p) where p = cos (46) 
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and the expression for dP/ (p)/af/U which is given in JA.45L we put Yf L in its final form: 

CO CO 

T n (r, 9) = ^ T, (r) sin 2 9 P, (cos ff) where T, (r) = J] rZr^k (r) 
£ [x] is the integer part of jc, and -Zjj.^ is defined in (1181 and spelled out in explicit form in Appendix B. 



4¥] 

2 (2fe-4j-l)%_ 2; _i 



(47) 



4.1. Evolution equation for the angular mean angular velocity 

Taking the horizontal average of equation J44> over an isobar and using the assumption that £2(r) » Q/(r), we obtain the following 
vertical advection/diffusion equation for the mean angular velocity CI: 

p^(r 2 Q) = ±-d r (pr 4 Uu 2 ) + -\d r (pv v r 4 d r £) + T> £ (r) , (48) 

J^r.%(r,0)sin0d0 



5r 2 



where the mean Lorentz torque is T<j= z (r) = 



sin 3 9A9 



r (r)--r 2 (r) 



(49) 



Note that only the / = 2 component of the circulation is able to advect a net amount of angular momentum; the higher order 
components of Um (for instance those induced in its tachocline by a differentially rotating convection zone) do not contribute to 
the vertical transport of angular momentum, as was pointed out in Spiegel and Zahn (1992). The explicit form of J48> . keeping 
only the I — {1, 2, 3} terms in the expansion of B (cf. eq. (II li ) is derived in Appendix E.2.1. . 



4.2. Evolution equation for the differential rotation in latitude 

We establish the equation governing the horizontal transport of angular momentum by multiplying eq. (I48> through sin 2 9 and 
subtracting it from the original form J44l >: 

p4 (r 2 sin 2 6>q) + V • (pr 2 sin 2 flflWj + ^-^-d r (pr 4 Q U 2 ) = ^V^<9,. (pv v r A d r Q) + -?—dg (pv h sin 3 9d e Ci) 
at v ; v ' 5r 2 v i r i \ i sin g \ / 

+ T^-sin 2 ^. (50) 

In the advection term we have again neglected the fluctuation Q(r, 9) compared to the mean £1 

The next step is to replace D(r, 9) by its expansion J32i in the horizontal functions Qi(9). For / = 2, the equation separates 
neatly into 

P j t (r 2 Q 2 ) - 2pQr [2V 2 - a(r)U 2 ] = ^d r {pv v r A d r Cl 2 ) - 10pv A Q 2 + T 2 (51) 



1 d , , v idln(r 2 Q) 

with V 2 = — — (pr 2 f/ 2 ) and a(f) = -—J ^. 

opr dr v ; 2 dine 



(52) 



It can be simplified by assuming that the turbulent transport is more efficient in the horizontal than in the vertical direction (i.e. 

p^ (r 2 Q 2 ) - 2pQr [2V 2 - orC/ 2 ] = -10pv,,Q 2 + T 2 . (53) 
In the asymptotic regime t » r 2 /v/,, a stationary state is reached: 

v h Q.2 = [2V 2 - at/ 2 ] Q + (54) 
5 lOp 

where horizontal diffusion balances horizontal advection and the torque due to the Lorentz force. 



For I > 2 the situation is more intricate, because there are couplings between terms of different /, which prevent a clean 
separation for each I. This is mainly due to the magnetic torque. Indeed, we recall that in the hydrodynamical case the hypothesis 
of Spiegel & Zahn (1992) allows such separation (cf. Paper I, §3.2.). Therefore, we choose here to stop the expansion of the 
rotation law at D 2 . The explicit form of J53i . keeping only the I = {1, 2, 3} terms in the expansion of B (cf. eq. (II li ) is derived in 
Appendix E.2.1. . 
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5. Structural properties of the differentially rotating magnetic star 

5.1. Baroclinic relation 

As in Paper I, we consider a non-uniform and a non-cylindrical rotation law, and add here a general axisymmetric magnetic 
field. Then neither the centrifugal force nor the Lorentz force derive from a potential, and therefore the isobars and the surfaces 
of constant density in general do not coincide. In order to determine how the density varies on an isobar, we start from the 
hydrostatic equation: 



1 Tlp 
-VP = g = -V</> + T C + — , 
P P 

where <p is the gravitational potential and g the local effective gravity, which includes both the centrifugal force 
T c = |^ 2 V (r 2 sin 2 o) and the meridional magnetic force T^p. Taking the curl of this equation, we get 

- —xVp A VP = — —Vp A g = iv(fi) 2 A V(r sin ff) 2 + V A |— P 
P P 2 \ p 

which to first order reduces to 



(55) 



P 



= [5 r (Q 2 )r cos 6» sin - <9 fl (Q 2 ) sin 2 fl]? f + VA 



P 



with p'(r, 6 1 ) being the variation of the density on the isobar. 

We now expand all terms in spherical harmonics. For the density fluctuation this is readily done: 



p (r,0 = £pKr)J»j(cos0). 



Then, differentiating the equation of state: 

dp dP dT dp 
— = a— -5— +<p— 
p P T p 

we can write the modal amplitude of the density fluctuation as 

BL = 0, = cpAi - 
P 

For the centrifugal force, we recall that 
Q(r, 6) = Q(r) + Y Q ; (r) (P,(cos 9) - /,) , 



which leads us to 



Q 2 (r,6») 



Q -2oYfii/, +2Q^O/P,(cos6») 

/>o J />0 

—2 . 

keeping only the terms linear in O; (beside £2 ). Finally, we draw the expansion of the poloidal Lorentz force from Jl 81 : 

00 

T £ ,p (r, 6) = Yj {Xr £ j(r)Pi(cos Gfir + ^r £ ;l(r)P] (cos ( 



(56) 



(57) 



(58) 



(59) 



(60) 



(61) 



(62) 



(63) 



It remains to insert these expansions in eq.|^l noting that it projects itself only on the azimuthal vectorial spherical harmonics 
(6>). Using the algebra related with these harmonics, we reach the following expression for the modal amplitudes of the relative 
fluctuation of density on an isobar: 



D,(r) + ± + -— r J 
rp(r) r ar \ p(r) 



(64) 



g is the horizontal average of the modulus of g and T>i (r) has been derived in Paper I (eq. 47): 
2>z(r) = Ni \rd r [if (r) - 20(r)0 2 (r)/ 2 -Ltf u 

+ 2rY J d r (Q(r)Q. s (r)) -1 [a° (-C°_A*-2 + D?_A S ) + (-C° +1 <5,, s + D° +1 <5 /)i+2 )] 

s>0 /Vj j 



(65) 
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where all the numerical coefficients (A?, Z?°, C?, D?, G°, i/, ) are given in Appendix A. 

This baroclinic equation J64l -I65> plays a key role in linking the density fluctuation on an isobar with the differential rotation 
and the magnetic field. It allows us to close the system formed by the induction equation, that for the transport of angular 
momentum, that for the transport of the chemical species and that for the transport of heat, which we shall establish in §6. 



5.2. Effective gravity and perturbed potential 



It remains to examine how the redistribution of mass we just described modifies the local gravity. In the absence of magnetic 

(66) 



field, the perturbing force is just the centrifugal force Tc - 5^ 2 V(r 2 sin 2 9), whose r and 6 components are expanded as 



r c M 8) = Yj a i(^ p i( cos r c,efr 0) = ~ J] bi(r)d g P,(cos 0). 

i i 

where a/(r) and bi(r) are given by: 

(<S/,o - 61,2) 



a,(r) = ^r[Q 2 (r)-2Q(r)Q 2 (r)/ 2 



+ 2rC2(r) V Q s (r) 



s>0 



b,(r) = ir[Q 2 (r)-2Q(r)Q 2 (r)/ 2 

+ 2rQ(r) J] ( A ° (-C^-Va + D^Afo.) + «! + ^ W )) • 

To this we have here to add the meridional components of the Lorentz force 
T £ , P ;r (r, 6) = J) As&i W Pj (cos 6) T l ,p; B = ~Yj M d o p i ( cos e ) • 



(67) 



(68) 



(69) 



In the expressions that we have derived in Paper I, it suffices then replace po fl / by po fl / + Xjzj ar, d Po^i by p^bi + to 
establish the version including the effect of the magnetic field. Thus for the gravity perturbation along an isobar we now get (cf. 
Paper I, eq. 72) 



gl 



dgo 1 



dr 



-zr\bi + 

So 



Pa 



+ — a/ + 

go 



Po 



d 

+ — 
dr 



\8o/ 



(70) 



The last term involves the fluctuation of the gravity field along the sphere, <pi, which is obtained by integrating the Poisson 
equation (cf. Paper I, eq. 60), modified along the same rules: 



1 d 2 / 1(1+1)- 4nGdp ~ 4/rG 

V* 1 ) - —^—^ -</» = — 



r dr 2 



go dr 



go 



d d / x 

poai + — (rpobi) + Xr £ ;l + j- \r^T L ;i) 



(71) 



Let us recall that Sweet (1950) was the first to establish this result for the most general perturbing force. 



6. Thermal imbalance and transport of heat 

It remains to implement the magnetic field in the heat equation: 



pT 



f +1/-VS 

dt 



V-(YVr)+pe- V-F h +J, 



(72) 



where S is the entropy per unit mass, x the thermal conductivity, e the nuclear energy production rate per unit mass and Fj, the 
flux carried in the horizontal direction by the anisotropic turbulence. Note that in a medium of varying composition, we have 
to take into account the entropy of mixing (cf. Maeder & Zahn 1998). In the simplest case, applicable to main-sequence stars, 
where the stellar material can be approximated by a mixture of hydrogen and helium with a fixed abundance of metals, it can be 
expressed in terms of the mean molecular weight only; then we have 



dS = C L 



dT dP dp 
TP p 



(73) 



where V a d is the adiabatic gradient and O is a function of the metal mass fraction and of p, the mean molecular weight. 
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The magnetic field manifests itself in the heat equation £721 through the Ohmic heating term J~, but also in the divergence of 
the thermal flux, because it involves the divergence of the perturbing force. This force includes both the centrifugal force Tq ar *d 
the Lorentz force per unit volume, iF^p/po, and their divergence will again be expanded in spherical harmonics as 



V-T c =fc + ^f c -> P,(cose) and W— )=/x + Z / ^ P ' (cos0) - 
;>o ' P° i i>o 

In Paper I (eq. 85) we gave the expressions of f c and fcj in terms of a\ and b\\ 
f c = ^<9,.(r 2 fl ) and f C j = ^3 r (r 2 a ; ) + 1(1+1)—. 



(74) 



(75) 



As above in 35.21 . it suffices to replace at by Xy^j/po and bi by J/rrj/Po to obtain the equivalent expressions for the divergence 
of the Lorentz force: 



and f £J = \d r I r 



,X 



+ 1(1 + 1) 



TC 



\ pq I r \ Po / rp 

This leads us to the modal form of the heat equation, which can be implemented directly into a stellar structure code: 
U,(r) 



<fPj ^dln/7 

+ <£ -A; 

df df 



■ Wad - V) 



M(r) p 



(76) 



(77) 



where T\ is given by: 

Jc + Jl { 



T, 



e + e„ 



1 - 

r 
3 







AnGp 



§l + foi±^i_fg±f^ ( _ sx¥i + (pAi) 



AnGp 



P 

(e+e. 



d r (H T d r ¥i -(1-6 + X TfV, - (f+X^i) 




AnGp 

1(1+ 1)H T 
3? 



{(H T d r % -(1-6 + X t) x ¥i - (f +^)A/) + (f e e T - fed + 6)% + (f 6 e, + f el p - V )A, 



(78) 



We recall that L is the luminosity, M the mass, T the horizontal average of the temperature, C p the specific heat at constant pres- 
sure, and V the radiative gradient. We have also introduced the temperature scale-height Hj = |dr/d In T\, the thermal diffusivity 

K = x/pCp, the horizontal eddy-diffusivity Dh and f e = e/(e + ? gra v), with e and e grav being respectively the mean nuclear and 
gravitational energy release rates, whereas and Xn WQ tne logaritmic derivatives of e and of the radiative conductivity x with 
respect to p., their derivatives with respect to T being noted as ej and^^ . Moreover, we have e m = L (r) /M (r) and p m is the mean 
density inside the considered level surface. 

Three remarks before we conclude this section. 

- First, we have noted previously that if the expansion of B is stopped at mode Z max = N m in M 1\ . then the expansions of Tz 
and of J' must include all terms up to Z max = 2N,„, due to selection rules. Therefore, keeping Q. (r, 6) — Q. (r) + Q.2 ( r ) Qi (#)> 
the expansion of the meridional flow will extend to ln M;J mx = max {4, 2Af„,}. 

- Next, the relative importance of the terms contributed respectively by the magnetic force and by the centrifugal force depends 
on the ratio (Va/Vq) 2 where Va = B/ yfp^po is the Alfven speed and Vq. = R£l, R being the radius of the star. If the magnetic 
field is weak enough, this ratio will be small except just below the surface because of the decay of the density. In this case, for 
the implementation in existing stellar evolution codes of the equations related to the thermal imbalance and to the structural 
properties of the star (cf. §5), we keep only the modes I = {2, 4), which allows to describe the tachocline(s) circulation to first 
order (see Appendix E). However, in the case of a strong field higher order modes of the Lorentz force must be taken into 
account. 

- Finally, the ratio between the term representing the Ohmic heating, J\ (r) /p, and the highest-order derivative term of V ■ (x VF) 
in j77!78> . namely (L/M)(p m /p)(r/3)d r (H T diyi'i), is given by (r//K) (Va/Vq) 2 - Therefore the Ohmic term can be neglected 
since Va < RQ and 77 <K K. 

In conclusion, the link between the circulation and its cause, namely the thermal imbalance due to the rotation, the magnetic 
field and the chemical composition, is established through J77l78i where the temperature fluctuation on an isobar is governed by 
an advection/diffusion equation from which we can derive the radial component of 14m- These equations have been established 
for a rotation law which depends on r and 9, and that allows us to treat simultaneously the bulk of the radiation zones and the 
tachoclines in presence of a general axisymmetric magnetic field. 
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7. Boundary conditions 

The system of equations is now complete: we have the induction equation, an advection/diffusion equation for the transport of 
angular momentum (mean and fluctuating) including the action of the magnetic torque, another for the temperature (mean and 
fluctuating), and the baroclinic relation which allows us to close the system. The equation for the transport of chemical species 
(or alternatively for the transport of molecular weight) is unchanged. It thus remains to specify the boundaries conditions of this 
system, to be applied on the limits of the radiation zone. To be specific, we consider a star with a radiation zone located between 
a convective core and an upper convection zone. We designate by r b and r, the respective radius of the base and of the top of that 
radiative zone. Of course, in a solar-type main-sequence star we have r b = 0, whereas for a massive main sequence star r, = R. 

The boundary conditions for the equations already present in the hydrodynamic case, which was examined in Paper I, are 
unchanged, except those to be applied to the equation of angular momentum transport, since it includes the Lorentz torque. The 
novelty here is the induction equation, which we have split in two equations, respectively for the poloidal and toroidal magnetic 
stream-functions & (eq. 1371 and %l ( e q- 1381 . These equations are of second order in r, and therefore each of them requires 
two boundary conditions. In principle they can be obtained by solving the induction equation, together with the equation of 
momentum, etc. in each adjacent convective zone, but this represents a formidable task which is well beyond our scope here. 
Instead, we make the following - simplifying but reasonable - hypothesis, valid only in the convection zone(s): 

- (i): we do not take into account the meridional velocity field. The only motion we retain is the zonal flow associated with the 
differential rotation: 11 ' v - r sin 6Q.cz (f, ff)^, with £lcz(r, 0) being the angular velocity in the convection zone; 

- (ii): we put rj v = r\h = ^turb where 77 tur b is the magnetic eddy-diffusivity, which we assume constant inside the considered 
convective region; 

- (iii): this eddy-diffusivity is high enough to allow for a stationary state (compared to the slow evolution in the radiation zone). 
With these assumptions, the poloidal field is potential, and d37i reduces to: 

it admits the general solution: 

4=Ar' +1 + ^. (80) 
r 

When applying it to a convective core where we must take the non-singular solution at the origin: 

r^-(l+ l)^o=0 at r = r b , (81) 
dr 

which becomes Q = if r b = 0. In the case of a convective envelope, surrounded by vacuum, we have 

r— + lA = at r = r t . (82) 
dr 

With the same assumptions, the toroidal field obeys the simplified form of ( 1381 : 

d*o _ £0+1) x i Si&cz) where 5 , (0cz) = _±.[^ ( o cz) + 5r(r j/ Ad; ,(Q cz ))]. (83) 
dr 1 r 1 77tui-b N 

<^Ad;/ (Qcz) and J^Ad;/ (Qcz) are derived using the previous multipolar solutions for (cf. eq.l80L expanding Qcz in spherical 
functions like Qcz (r, 0) = 2/ ( r ) Pi ( cos where the £2/ are taken from the results of the inversion of helioseismological data 
(Corbard et al. 2002), using J34I and Appendix C. Then, the solution of fl83i inside a convection zone is derived using Green's 
functions: 



1 1 1 



Xo = 



21 + 1 77turb 



1 x- l Si(Qez) djc + r-' x M S, (Qcz) dx 

Jr •JRmlCZ 



(84) 



where Rj„,cz and R ext cz are respectively the radius at its base and at its top. Applying (l84i to a convective core where /?,„,cz = 
and flextcz = r b , we get: 

Xo = „/ . 1 r~ b l [ " x l+1 S, (Q cz ) dx at r = r h . (85) 

11+ 1 7/tui-b J() 

Finally, for a convective envelope where R m ,cz = n and R ex tcz = R, we obtain: 

^o = ^T — r \ +l f *"'S/(Q CZ )d* at r = r t . (86) 
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In the case where Si (£2cz) is negligible, the previous results leads to: 
Xq = at r — ri, and r — r,, and thus to Bj = 0. 



(87) 



It remains to write down the boundary conditions to be applied on the transport of angular momentum. Since the equation for 
the mean angular velocity Q. d48l is of second order in r, it requires two boundary conditions. They are obtained by evaluating 
the budget of angular momentum in each adjacent convective zone: 



d 
df 



r 

Jo 



r 4 pQdr 



1 a — d 

-rpQ.U 2 - T B (r b ) at r = r b , — 
5 df 



r 



r 4 pQdr 



l -r 4 pW 2 -Ta + r B {r t ) ntr=r t . 



(88) 



Tn is the (signed) flux of angular momentum which is lost at the surface by the stellar wind, and Tsir) the flux carried by the 
magnetic field through the considered surface: 



r B (r) = --r 3 f B r (r, 6) (r, 0) sin 2 GdO = -L V \l (I + 1) JV% 
HO Jo ^tol 



(21+ l)(2/ + 3) Xo 



(2i- i)(2/+ir° 



(89) 



Here again the perturbation Q. 2 obeys an evolution equation which does not include any derivative in r, and therefore it needs no 
boundary condition. 



8. Conclusion 

The work presented here is the continuation of Paper I (Mathis & Zahn 2004), where the star was assumed without magnetic 
field. Here we allow for a magnetic field of moderate strength in the radiation zone(s), with an Alfven speed not exceeding the 
rotational velocity: (Va) 2 < (R£l) 2 . Such a field has little impact on the hydrostatic balance, since we also assume, as in Paper I, 
that the centrifugal force is small compared to gravity: RQ. 2 <K g. But it will compete with the centrifugal force in the baroclinic 
balance d56l >. and therefore it will participate in governing the meridional flow. Moreover, such a field tends to enforce uniform 
rotation along the field lines of the poloidal field (Ferraro's law). The question then arises whether or not the poloidal field threads 
into the convection zone(s), where such differential rotation is probably maintained through the turbulent motions, as observed 
in the Sun. To answer that question one has to treat consistently the evolution of both angular velocity and magnetic field, with 
special care for the dynamics in the tachocline(s), and for this reason we expanded them in spherical harmonics to arbitrary order. 

The overall problem of rotational mixing in a magnetized star is highly non-linear, with multiple feed-backs, as illustrated 
by the diagram displayed in fig-E Mixing is achieved through the meridional circulation, which is due to the thermal imbalance 




Fig. 1. Rotational mixing in stellar interiors with magnetic field: a highly non-linear problem (the action of the convection through 
penetration and overshoot which enhance the mixing has been added). 



14 



S. Mathis and J.-P. Zahn: Transport and mixing in rotating stars: II-Axisymmetric field 



caused by the centrifugal force and by the Lorentz force (mainly through the baroclinic balance), and also through the turbulence 
generated by the shear of differential rotation and possibly by the magnetic field. These motions modify both the rotation profile 
and the magnetic field through large-scale advection and turbulent diffusion. As the star evolves, molecular gradients build up 
through microscopic diffusion (including radiative levitation and gravitational settling), and through nuclear burning. Large-scale 
advection turns them into horizontal gradients, which react on the meridional circulation and are smoothed by the turbulence. 

The boundary conditions play a crucial role: in a non-magnetic star, the circulation is driven by the loss (or gain) of angular 
momentum, as explained in Zahn (1992); in the absence of such loss, the circulation tends to vanish, as was shown by Busse 
(1982). It remains to be seen whether this is still the case when angular momentum is transported mainly through magnetic 
stresses, with the poloidal field anchored in the convection zone(s). 

The main weakness of our modeling remains the description of the turbulence. As in Paper I, it is assumed to be anisotropic, 
due to the stable stratification, and that it tends to smooth angular velocity and chemical composition on horizontal surfaces. The 
action of magnetic field on such turbulence is not taken in account, which seems reasonable as long as the Alfven speed does not 
exceed the rotational velocity. Moreover, we have deliberately ignored here the possibility of MHD instabilities generating their 
own turbulence (Balbus & Hawley 1994; Spruit 1999). 

From the technical point of view, we adopted here the formalism developed in Paper I and introduced the spherical vectorial 
harmonics, which obey the Racah-Wigner angular momentum algebra widely used in quantum mechanics. This allowed us to 
separate the colatitude 6 in the vectorial partial differential equations which govern the problem, and to reduce the problem to 
solve coupled partial differential equations in t and r only. It is then straightforward to introduce these equations in existing stellar 
structure codes, a task we have now undertaken. 

Among the first applications we plan to model a solar-type star, to check whether a fossil magnetic field, which is able to 
extract angular momentum from the radiative interior as the star spins down, yields a profile of angular velocity that is compatible 
with the helioseismic data. Next we will undertake the modeling of magnetized massive stars, since they are progenitors of 
neutron stars which are the seat of strong magnetic fields (Maeder & Meynet 2003, 2004; Heger et al. 2004). 

Most of the Appendix is dedicated to the algebra involving the spherical and vectorial spherical harmonics used in this paper. 
The numerical values of the coupling coefficients given in the Tables have been computed with Mathematica. In Appendix E, we 
state the equations in a form ready to be implemented in a stellar structure code, when the magnetic field is expanded up to the 
octupole (/ = 3) and only the first term D.2 of the differential rotation is retained. 

Acknowledgements. This work was supported by the Centre National de la Recherche Scientifique (Programme National de Physique Stellaire 
and GdR Dynamo). 
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Appendix A: Algebra related to the spherical harmonics 

A .1. Scalar quantities 

A.1 .1 . Definition and basic properties 

The spherical harmonics are defined by: 

Y',"(0, <p) = NfP^ (cos 0) e"^, (A. 1) 

where f| m ' (cos 0) e""^ is the associated Legendre function, and the normalization coefficient being 
21+ 1 (Z-|m|)f ' 



(m+H) 

= (-!)— 



4-n (I + \m\) 



(A.2) 



They obey the orthogonality relation: 

£ (y;;" (e, <pj)' (e, & do. = 5 hh 5 mum2 (a.3) 

where dQ = smOdOdip and where the complex conjugate spherical harmonic is given by: 

(r?(fi,<pjj = (-\) m Yj m {e,<p). (A.4) 

Using these properties, every function f(6, <p) can be expanded as: 

/ (0, <P) = 2 2 f!,Xi & <f) where f m = f(0, <p) (y? 1 (6, <pj)' d£2. (A.5) 

1=0 m=-l J Q 

A.1.2. Special case of axisymmetric spherical harmonics 

In that case we have: 

If (0, if) = NfP l (cos 0) and d e Y? (0, tp) = N°d e P t (cos 0) = j ° j^i = ( £ s g) for l > Q (A.6) 

A.1.3. Linear differential and recursion relations 

We recall that the spherical harmonics obey to the differential equation: 

—d e (sin 0d e Y? (0, ip)) + — (6. V) = -/(/+ 1) If ^ V) ■ (A.7) 
sin0 v ' sin 2 6» v 

In deriving the functions which are related to the centrifugal force in §5 (cf. J65I67I681 ) we have used the following recursion 



relations for m = 0: 



cos 0Yf(0) = AW , (0) + BW.A6) where A? = 1 — and B° = . (A.8) 

' ' ' 1 ' /+1 ' V(2/+ 1)(2Z- 1) V(2/ + 3)(2/+ 1) 
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sin 9Y?(9) = C^dgY? , (9) - D°d g Y° (9) where C? = 1 — and D? = 1 . (A.9) 

' 111 11+1 ' V(2/+l)(2/-l) VC2T+3X27TT) 

c„ s *y,» ((! , = where £l U _^L_ a „ d ^ = _L_, (A.,0) 

-of*, «,-!*»•,«) where of. and H? - — g ^ _ ■ ■ . (A.H, 
The identities ( IA.9> and (|^^) have been deduced from the two others jA.8IA.lH with the help of 1A.7I . 
A.1 .4. Expansion of products of the spherical harmonics 

Using the normalization and the orthogonality of spherical harmonics (cf. jA.3IA.5t ) and their complex conjugate (cf. jA.4» . we 
can write: 



Y"' 1 (9, <p) Y™ 2 (9, ip) = (-l)(»»+»^ ^ Z™ 1 j%'~ < - mi+m * ) Yj' 1+m2 (9, <p) (A.12) 

i=\h-h\ 

where we define the integral J™',;7' m like in Edmonds (1968) or Varshalovich and al. (1975): 

isr- X»r«rtT«rtT«rt->- V ' 2 '''"' 2 '::'" 2 " " (,:.', t J)(5 o o) < A13 > 

with the 3j-Wigner coefficients which are related to the classical Clebsch-Gordan coefficients by: 

\m, m 2 mj y 2 / + 1 l i,m\,h,m 2 - v ' ; 

So finally, we get the following expansion for the product of two spherical harmonics: 

h+k 

Y?(0,<p)Y£{0,<p)= J] c' Mm Y?^(e,<p) (A.15) 
i=\h-h\ 

where the coefficient c\ , is given by: 



J -r_n(m l+m2) l(2h + l)(2l 2 + l)(2l+l)( 1, h I \(hhl 



c W/ 2 ,»>2 ( 1)m ' '" 2 "V 47r (m, m 2 -(m, +m 2 ) j( j' (A ' 16) 

Then, using the initial definition of spherical harmonics (cf. IA.HA.2t . we deduce the expansion for the product of two associated 
Legendre functions: 

P'/;' (cos 9) (cos 9) = £ <, m „ /2 , m ^7 1+m2 (cos 0) (A.17) 

/=l'i-fcl 

where 



d l =(-U^ + »»U2l + D (h+m 1 y.(l 2+ m 2 y.(l-(m l+ m 2 ))! l h l 2 I \(hl 2 l\ 

%, mu i 2 ,m 2 k v + -y ft _ mi) , (fc _ otz)! (/ + (mi + m) y \ m m2 _ (OTl + m 2 ) j ^ o o o j ■ lAa6J 

Note that the previous expression is symmetric; therefore 

^\\,m\h,mi ~ ^\ 2 ,m 2 ,h,m\' (A. 19) 
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A.2. Vector fields 

A. 2.1 . Definitions and basic properties 

Following Rieutord (1987), we expand any vector field u(r, 8, <p) in vectorial spherical harmonics as 
u(r, 8, <f>) = J] J] {u l m (r)R™(G, <p) + v l m (r)S?(8, <p) + w l m (f)T™(6, ?)} , 

1=0 m=-l 

where the vectorial spherical harmonics R"' (8, <p), Sf (8, <p), T"' (8, <p) are defined as: 
Rf(8, cp) = Y™{8, <pje r , Sf (8, <p) = V S Y?{6, <p) and I* (8, <p) — V$ A R?(8, <p), 

with the horizontal gradient = Tgdg +T IP cL . 

sm6 

These vector functions obey the following orthogonality relations: 

f R™ 1 -S^dQ. = f R? ! ■ T"' 2 dQ. = f S? 1 • Tfda = 0, 
Jn 1 2 Jn 1 2 Jn 1 2 

Jf jjm, . = $ ljfe <V,m 2 and £5;;' ■ (s'ffdO. = f n K- { T Z)* dn = h(h + l)S h ,hS mum2 . 

The vector function u may also be projected on the classical spherical vectorial basis: 
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OO I 



Z=0 m=-l 



e g + 



v ™ (r) ^W v) - wL( <p) 

sin a 



(A.20) 
(A.21) 

(A.22) 

(A.23) 
(A.24) 

TJ, (A.25) 



T r , Tg and^ beeing the unit- vectors respectively in the r, 8 and tp directions. These expansions of vector fields allow us to separate 
explicitly the angular variables 8 and tp in the vectorial partial differential equations which govern the problem. We thus reduce 
the problem to solve partial differential equations in t and r only, which are easy to implement in existing stellar structure codes. 
Note that Rf {8, <p) and 5™ (8, <p) represent the poloidal part, and T'" (8, <p) the toroidal part of u. 



A. 2. 2. Expansions of differential operators 

As stated in the previous section, the expansion of the vector fields in Rf (8, <p), S"' (8, <p) and T" (8, tp) allows us to separate 
the variables in the vectorial partial differential equations which govern the problem. This will prove particularly useful when 
dealing with magnetic field, and with the non-linear expressions where it is involved. We start by expanding the classical linear 
vectorial operators: gradient, divergence, curl and laplacian (scalar or vectorial). 



Gradient: 

Taking a scalar function f(r, 8, <p) expanded in the Y" 1 (8, ip): 

CO / 

1=0 m=-l 

it is straightforward to derive its gradient: 
V/(r, 9,0 = ^ d r &R!?(8, tp) + J -^Sf{8, tp) 

1=0 m=-l ^ r 

Divergence: 



(A.26) 



(A.27) 



Taking a vector field u expanded as in (IA.20> . its divergence is given by: 

l r , J 



1 V 

V-u(r,8,<p) = Y,T; -dr{r 2 u' m )- 1(1+1)^ 



1=0 m=-l 



Y?(8,<p), 



(A.28) 



where whe note that, if u is divergence-free, such as the magnetic field B or the momentum density pliu m the anelastic 
approximation, we have the following relation between u' m and v' m : 



v — 



(A.29) 
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Laplacian of a scalar function: 

With these expressions for the gradient of a scalar quantity and for the divergence of a vector field, one can easily derive the 
laplacian of a scalar function. Using the well-known property V 2 / = V ■ (V/), we get: 



v 2 / = v • (v/) = J] J] -d r > (rf m ) - i(i + i)f r;\e, = ^^ Aj^e, V ) 



1=0 m=-l 

where A/ is the laplacian operator: 

A; = 8 r ,r+ -8 r 5 ■ 



/=0 m=-l 



(A.30) 



(A.31) 



Curl: 

This operator is found in the expression of the magnetic field in terms of a stream function, in that of the current density, etc. 
If we take u expanded as in ( IA.20t >. we retrieve Rieutord's (1987) result: 



1(1+1)- 



R?(0,<p) 



d r (rw\ n ) 



sT(o,<p) 



uL 1 



dr(rv l m ) 



r;'(0,<p) 



1=0 m=-l 

Taking this result, the expansion for the laplacian of a vector field is derived. First, we have: 

OO / 

1 I II I +11/ II I , , \ 



V A (V A u) = J] 



1=0 m=-l 



*->K)) 



-d r u' m - -d v (rv' m ) 
r r 



Sf(0,0 + [-AyJrf (6,<p)\ (A.33) 



and therefore 



V 2 m = V(V-h)- V A(VAh) 



°° ' IT 2 

A '"m-^K-/a+DvL) 

/=() ;«=-/ " ' 



R"'(0,<p) + 



Aiv m +2~ T 



ST(0,ip) + [A,w l m ]Tf (6,^)1, 



result which becomes in the case where u is divergence-free (cf. JA.29> ): 



V 2 H = -V A (V A u) 



IE 



1 



-a,K,) 



1 / A,(r<) 
r T /(7 + 1) 



(=0 m=-i 

Finally, using jA.32> once again, we get the 'triple curl' operator: 



Sf(6, l p) + [A l w' m ]Tf(e,^ 



oo / 



1=0 m=-l 



-1(1 + 1) 



A,< 



--d r (rA,*i) S?(8, <p) + [A lZ ' m ] 77(0, 



where: 



;5r(rvL) 

which becomes in the case where u is divergence-free (cf. < IA.29I > ): 



CO / f 

(VA^^W 

1=0 m=-l 



-1(1+1) 



A,w, 



R?(6,<p) 



--d r (rAiw' m ) 



S?(8,<p) 



1 



1(1 + 1) 



AiAi(ru' m ) 



T?(6,<p) 



(A.32) 



(A.34) 



(A.35) 



(A.36) 



(A.37) 



(A.38) 



A. 2. 3. Products of axisymmetric vectorial spherical harmonics 

Before we turn to the non-linear terms involving the magnetic field: the advection term in the induction equation and the Lorentz 
force, we shall first perform the projections of the scalar product of two general axisymmetric vectors and of their vector product, 
and derive the associated coupling coefficients. 



Scalar product: 
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We take two general axisymmetric vectors X\ (r, 9) and X 2 (r, 9), which we expand as in JA.20> : 



Xi (r, 0) = £~ =0 b§ (r)*° (9) + fig (r)S° (0) + Cj (r)J° (0)} 
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(A.39) 



and perform their scalar product. After some algebra involving the expansion of products of spherical harmonics (cf. §A.1.4.), 
we obtain 



X x (r, 9) ■ X 2 (r, 6) = £ P (Jf , . ZlW (r) P, (cos 0) 



with the following expression for P(x l x 2 );i( r )'- 



CO CO 



h+k 



I j=/(; l; o,/ 2; o) 



/,=o/,=o 



(A.40) 



(A.41) 



Here ; is defined in terms of the coupling coefficients d\ m l m (IA.18> : 



/,+/, 



X 



- Z 



09 



<i,« Z [<2./-4/>-l)-X] 



p=0 



(A.42) 



with 



2 [2(./-2/,- 1)- 4^-1] 



9=0 



[a-2p-l)-2 ? -l]+2 

^.[O^-D^-l] - Z 

r=/(2,0,[O-2/7-l)-2 ? -l],0) 



f 2,0,[y-2p-l)-2 ? -l],0 (5 '. r 



(A.43) 



We have used the classical notation £y for the usual Kronecker symbol, E [x] is the integer part of x and 
I(l\,ni\, I2, mi) — max — l 2 \, mi + m2). 

Vector product: 

We operate likewise for the vector product of two general axisymmetric vectors X\ (r, 9) and X2 (r, 9), again expanded as in 
JA.20> . We reach the following result: 

CO 

Xi(r, 9) A X 2 (r, 9) = X iX ^x 2 ydr)? r + J] [X( Xl hX^{r)Pi (cos 9)7 r + V ' ( x^x 2 yi(r)P\ (cos 9)7 8 + Z(. Xl AX 2 y,i{r)P\ (cos 0)7^) 



/=() 



X(x,AX 2 y,i( r ) 



-^(X 1 AXi);l( r ) 



Z(x l AX 2 yA r ) 



where we have also used the following property of the Legendre function (Gradshteyn & Ryzhik 1967): 



dp 



\ z / 

J] (2l-4k-l)P l .2k-M- 



The radial function X(x X AX 2 yj (r) is given by: 



(A.44) 



(A.45) 



(A.46) 



/i=0/i=0 



with the same coefficient Xl , derived above in JA.42I I. The two other radial functions if(x,AX 2 yi ( r ) an d Z^ca'^a:,) / ( r ) are given 
by 



CO co 



/,+/: 



J/ (M); /(r) = J] J] [c*(r)£#(r)] £ <ia,<Aj " MW^(r)] Z <0A.i tf « 



and 



/i=0/ 2 =0 



00 do 



i=/ft,o,; 2> i) 



(A.47) 



(A.48) 



/i=o; 2 =o 



;=/(/,, l,/ 2 ,0) 



J=I(h,0,k,l) 
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Appendix B: Lorentz force 

In this section we shall project the Lorentz force, 



-(VAB) 



AB, 



(B.l) 



on spherical vectorial harmonics, and write the result in terms of the poloidal and the toroidal magnetic stream-functions, respec- 
tively and^-g. We have just seen in §A.2.3. how to expand the vector product of two axisymmetric vectors in the R l J(0), 5^(0), 
T°(0). We apply the method to flTTt . with (cf. (O, O) 



1 °° 

Xi (r, 6) = j (r, 8) = -(VAB (r, 6)) = £ (r)fl° (0) + fig (r)5° (0) + eg (r)J° (0)} 



1 00 



Zi=l 



r9(fl) 



and (cf. O) 

DO 

Z 2 (r, 6) - B (r, 0) = £ (£#(r)*° (0) + S{, 2 (r)S° (0) + ^(r)?^)} 



Z 

fa=i 



+ 1)4 



5?(0) + 



fa 
AT 



So, we have, using the definition of §A.2.3. 



(B.2) 



(B.3) 





r>'i _ i 



-A, 



and 



T) 2 = 



fettz + l)f 



4 2 - 



~ rt) 

Therefore, we get using (IA.44IA.46IA.47IA.4gt : 

7£(r,0) =Xi(r,e)AX 2 (r,6) 

oo 

= fe>(r)e, + J] (fe(r)J», (cos 07, + feWP, 1 ( cos ^ + feW^/ 1 (cos 0)Tj 



/=! 



A' 



7\c;0 



A? 



*?($)• 



-few 



s?(0) + 



feW 



T l !(0)\, 



(B.4) 



(B.5) 



where the radial functions Xj^ j (r), (r) and -Z^;/ W m ' e: 



2 1 



OO oo 



3yU 



^4 



fa 



/i./ 2 



(B.6) 



*«=^zz>« 



;,=t fa=t 



-A,, 



^ eh \ 
?0 



A 

Uh + 1)4 



\ h+h 

X d hXh,0 S U ~ 

j=I(h, l,fa,0) 



r 



Zi+fa 



Z <0A.lM (B ' 7) 
./=/(/, ,o,/ 2 ,n 



CO CO 



^ u ;. = i ;„=1 U 



wfc + 1)4 



\ 'i+fa 



j=/(/i,0,fa,D 



Explicit values for X\ l if^^E) are given in Table lB.il and in Table IB. 2l for the coupling coefficients d\ Y ; . 
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Table B.l. Values of the X\ l coefficient involved in the r-component of the Lorentz force and in the Ohmic heating 



X h,h 



h \ h 


1 


2 


3 


h \ k 


1 


2 


3 


1 


-1 





26 
7 


1 








18 
35 


2 





93 
35 





2 





108 
175 





3 


2d 

7 





478 
49 


3 


18 
35 





28898 
2695 



Table B.2. Values of the d 1 ^ x ; Q coefficients involved in the 0— and in the ip— components of the Lorentz force and in the Ohmic 
heating 







(1 




d 2 
h,l,h,0 


d 3 
h, 1,(2,0 


h\h 


1 


2 


3 


h\h 


1 


2 


3 




1 


2 


3 


1 





i 

5 





1 


i 

3 





i 

7 


1 





i 

5 





2 


3 
5 





9 
35 


2 





1 

7 





2 


2 
5 





i 

15 


3 





18 

35 





3 


4 
7 





2 
21 


3 





1 

5 





d A 


d 5 

/],l,/ 2 ,0 




h\h 


1 


2 


3 


h \ h 


1 


2 


3 




1 


2 


3 


1 








i 

7 


1 











1 











2 





9 
35 





2 








4 
21 


2 











3 


3 
7 





9 

77 


3 





2 

7 





3 








50 
231 



Appendix C: Advection term of the induction equation 

We deal likewise with the advection term in the induction equation: 

(I/^+I/mJaB, (C.l) 
where ll^r, 6) = r sin 6Q.(r, OJF^, and Km represents the meridional flow, whose expansions are found in d30l . Therefore 

DO 

X l (r, 6) = tl v (r, 6) + U M (r, ff) = J] (0) + (r)S° (0) + C* (r)2* (0)} 

CO 

= 2 K 1 m < ( ) + v o w s ° ( ) + w 'd w r ° 



i,=0 



(C2) 
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Table C.l. Values of the X| ^ coefficients which are involved in the r-component of the induction equation 



X h,h 



X lh 



1 

2 
3 
4 



^ 



x2 

35 



13 n 366 
U 49 



3594 
245 



/, \ h 
1 

2 
3 
4 



^ 



26 

7 



26 n 478 



226 
21 



/l \/ 2 
1 

2 
3 
4 



f 



22 

15 



o s 





774082 
38115 



and (cf. (QT)) 

oo 

1 



'2=1 



-s- 



r 2 



r 



S?(0) + 



So, we have, expliciting the velocity field: 



^0 "0 JV° 

S'' = v' 1 = ^ 

■°0 w° 



Cg = wg = r 



1 c 4 , 
— o* Hio* 



. O* (r) = O (r) + ifi 2 (r) 
where { °; ; _ , \ y . 5 , n and ^ 
Qt (r) = Q/ (r) for / > 



Hence we get, using IA.44IA.46IA.47IA.4fl : 

(l/„(r, 0) + l/j^Cr, 0)) A B(r, 0) = Xi(r, 9) A X 2 (r, 6>) 

00 

{^Ad;l(r)P/ (COS 0) ? r + J/ A d;l(r)P, 1 (COS 0)^ + ZaA^P] (COS 6)7 v } 



1=1 



z 

1=] 



*Ad;l(r) 



-#Ad;l(r) 



A? 



where <Y A d ; i (/), J/Ad ; i W and ZAd ; i 0): 



ft=0 i 2 =l 



-^4 



z 2 



X 



D 2 = 



4 2 = \:M] 



J 



CO CO 



ft =0fa = l 



ft +fc 



X d hXh,0 6l >j 

) j=I(h,U2,0) 



ft -^0_ 

r 



h+h 



Yj d Lo,h,i 5l <J\ 



CO CO 



Zw = YY« 

'i =0/2=1 



h(h + 1) 



h+h 



> j=l{h,\,h,G) V 



Wo 



>=/(/,,0,/,,l) I 



Explicit values for X| l (|^^} are given in Table KJ.fl and in Tables KTJ2"llC.3IKT4l for the coupling coefficients 
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Table C.2. Values of the coupling coefficients which are involved in the 9- and in the <p- components of the induction equation 
for I = 1 



d' 

h.Uh.O 


d 1 

/i ,oy 2 .i 


h \ h 1 2 3 

1 -\ 

2 a o 

Z 5 " 35 

3 if 

4 f 


h\h 1 2 3 

1 f 

2 -i iS 

A 5 " 35 

3 0-^0 

4 -f 



Table C.3. Values of the coefficients involved in the 9- and in the tp- components of the induction equation for 1 = 2 



,2 


d 2 

il ,0,'2,1 


h\l 2 1 2 3 

1 i -i 

2 i 

3 i i 

7 " 21 

4 0^0 


h\l 2 1 2 3 

1 I i 

3 " 7 

2 i 

3 -7 2T 

4 -A 



Table C.4. Same as in TableOl but for I = 3 



d 3 

'l,l,/2,0 


d 3 


Zi \ / 2 1 2 3 

1 i 

2 2 o -i 

z 5 u 15 

3 | 

4 - — 

"* 9 " 33 


h\k 1 2 3 

1 f 

2 i 1 

z 5 " 5 

3 0-^0 

4 -i i 

9 " 33 



Appendix D: Ohmic heating 

For sake of completeness, we also calculate the Ohmic heating rate, which we shall express in terms of the magnetic stream 
functions A and^[ } . From §2.4, we get: 

J=— [IMI®(VA£)]-(VA£), (D.l) 
/Jo 
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where we allow for different eddy-diffusivities (77/,, 77,,) respectively in the horizontal and vertical direction. We apply again the 
method of §A.2.3., and identify the two vector functions which enter the scalar product: 



1 DO 

X l (r, 6) = - [\\tj\\ ® (V A B(r, 6))] = £ (r)«° (0) + flg (r) S° (0) + C§ (r) J° (0)} 

r9(fl) 



'i=i 



1 y"0 



r 



1 /A) 



-A,, 



( eh \ 
So 



and 



Z 2 (r, 0) = V A B (r, 0) = £ (r)«° (0) + 4 2 (r) S° (0) + (r) J° (0)} 



-S- 

z 2 =i 



wz2+i)4 



< 2 w + 



r 



Si (0) + 



( eh \ 
So 



We thus obtain, using JA.40IA.4T1 i: 

CO 

S (r, 0) = Z, (r, 0) • Z 2 (r, 0) = J] Ji to p > ( cos ) 



with the following expression for Si ( r ) : 



U=l fa=l 



2^ 
3 no 



-Mo 



h(h + A 



( eh \ 

So 



k 

h(h + 1)^- 



\ h+h 



So 
r 



x: 



luh 



(D.2) 



(D.3) 



(D.4) 



(D.5) 



recalling that X| ; is given in IA.42L and in explicit form in Table IB. fl 

Appendix E: Equations to be implemented in stellar evolution codes 

In this section we shall give the equations ready to be implemented in a stellar structure code, with all coupling coefficients being 
replaced by their explicit value, in the special case where only the dipole (I - 1), quadrupole (I = 2) and octupole (I = 3) are kept 
in the magnetic field, and where the differental rotation is reduced to its first term (1=2). 

E.1. Induction equation 

The linear terms of this equation readily separate into their poloidal and toroidal components, and they project on a single 
multipole, whereas the advection term is the result of various couplings. 



E.1 .1 . Equations for the dipole 

The induction equation translates into two evolution equations for the magnetic stream functions: 
^ = 2 ^| rZ Ad;I + ifrrAi (f ) 



it - - x ; i'-<-Ad;l + Tjh 1 



d/ , „ r v //t „ - [*Ad;l + r (rJ/ A d;l)] + 

and where the advective terms are given by: 



2 2 
where Ai = d nr H — d r — r, 



' = j± ji(cisi - Bp*) + f ^(cisi - vm) - h (123 ^Kn + 5 -f Kn) 



^Ad;l - 5p 



^Ad;l - JZ 



Vf ( 3S o^o + - -k ( S o»o + ^o) + Ti ( 5 W> + 3 ^o fi o) 



(E.l) 



(E.2) 
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with 

3§ = \ V5FC/4 
S 2 n = 2^V 2 



c 





Ml 


£> 2 = 


6% 






So = [754] and • 




T 1 = 

7 




T 2 = 
7 


Xo 
r 





" £ 3 " 



E.1 .2. Equations for the quadrupole 

Likewise 

f =2^/frZ A(1 ;2+^A 2 (|) 

and 

<^Ad;2 



2 [*Ad;2 + dr O^Ad*)] + <?r (l/A^) - 



where A2 = 3 r , r H — <9 r 



2 j_ 

3 2* 



#Ad;2 - 5p 
-3Ad;2 = 7= 



13 (C„<5 ( , + Cq£q) - §CjfiJ + ^Cl&l - TsBffij; - ™ ^K^o 
f (s 2 © 2 - y[ ( 2 fi 2 ) + f V5 (SflJCg + 2^£ 2 )] . 



All functions Jft? . . . Tq have the same meaning as in jE.3> . 



E.1 .3. Equations for the octupole 

The result is similar for the octupole: 

'f =2VfrZAd ; 3+wA 3 (f) 

^ + <9 r (r)^ = 2 Vf [*Ad;3 + <? r (ri/ A d;3)] + \d r (mdrXo) ~ 



2 12 

where A3 = d rr H — 3 r 7- 



387041 1 



and 

^Ad;3 = § j [| Vf (C, 1 ^ 2 + Sg^ 1 ) + 12 (C ( 3 )£ 2 - S 2 ^ 3 ) - f V^S^ 1 - ? _^ n , n 

^Ad;3 = ± I Vf " ^) + VH C 0^0 " ^) + ^0%' " TT^O^ 



z 



Ad;3 - 4 „ 



Vf ( 2 ^o - + ^ (5^1 + + 4 " ^o) + TT " ^0 £ 



E.2. Mean equations 
E.2.1. Mean rotation rate 

We recall the transport equation for the mean angular momentum J48l >: 

where the mean magnetic torque is given by: 

Yr £ = To - -T 2 with \ v , v ' 

with 

' = i { Vf [3 (^o " ^o) " " + if [2 (^o - .flgSg) - (figflg - ^fig)]} 

• ^ = 47 { Vf [ 2 - + {K&o - + Vf [( S o»o - ^o) + 5 (^o - ^o fi o)]} 

. Znfi = 27 V^ [( S o^o - ^o) + I ( S o»o - ^o £ o)] ■ 
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The St, &L C' DL, 6' and T' with / = {1, 2, 3} are given by: 



<7fl - 1 

r 1 - J- 



0' w 

24 



-Ai — 



*3 



2% 



<z,2 _ 1 



/Jo 

C 2 - ±- 



6^ 



-A 2 a 



f 2 
6§ 



[>4l 



and 



^ 3 



C 3 = J- 



From ( IE.9> . the final form of 1E.8I is immediately derived: 
Pj(r 2 Q.) = -L5,(pr 4 nt/ 2 ) + \d r (pv v r 4 d r Tl) + rZr £ ;i- 



f 3 

&4 



(E.ll) 



(E.12) 



E.2.2. Mean chemical composition 

We restate the mean transport equation for the concentration of chemical species, which remains unchanged by the introduction 
of the magnetic field (Paper I, eq. B.2): 



(E.13) 



E.3. System fori = 2 

E.3.1. Meridional circulation 

We split J77I78I in two first order equations as: 



U 2 = 



Mg 



\PC P T , 



1 



' ad 



(E.14) 



where 
S 2 = 2 



. fc + fz ( 
4nGp 



e + e 



grav 



# + fca + ft 2 - f -^U (~^2 + ^A 2 ) + ^ 
g 47rGp 47rGp p 



r „ 2H 
~d r ^2 - 
3 r 



*?5 



+ {^2 + if e e T - fed + 6)V 2 + (U, + f e <p - V )A 2 } + - 



J 2 -/d*F 2 dlnp \ 



and 



m 2 = H T d r ^ 2 - (i - 6 + XT m - Qp+Xr)M. 

The coefficients related respectively to the centrifugal force and to the Lorentz force are given by: 
7c = -Z^ao) , Tea = id, (r 2 a 2 ) + 6^, 



where 



PO 



and /x,2 = \d r I r 



2 ^?i;2 



Po 



+ 6- 



rpo 



a 2 = -|rO + |irOQ 2 
£ 2 = irQ 2 + ^rQQ 2 



(E.15) 
(E.16) 

(E.17) 
(E.18) 

(E.19) 



and 



= ~ i^f [(Cgfij - fifo 1 ) + (cjfig - S^)] + J (C £j - Sjrj) + | (cgfig - flgf*) + I (c 3 £ 3 - S 3 ^ 3 )} (E.20) 



x. 



r z ;2 - 



2 

3 2n 
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(E.21) 



^ = { Vf ^ ^ " ~ ^ ~ + ( C(1, ° ( » " ^ ) + 7 ^ " + \ ^ ~ ^ )f - (E ' 22) 



The C* Q , D l , &[ and :F ' with Z = {1,2,3} are given in (lETfY 

The relative fluctuation of the effective gravity is given in J70l >: we apply it here to I — 2: 



82 
g 



5 r \b 2 + 



— \a 2 + 



d 

+ d7 



Po / go \ Po 
where 02 is solution of the Poisson equation (cf. eq.l71>: 



(E.23) 



1 d 2 / - v 6 - 
r dr 2 v ' r l go dr 



4nG dpo— A-nG 
go 



Po«2 + -j- (rpob 2 ) + X.<f l -2 + — (rj/r £ a) 



dr 



(E.24) 



The Ohmic heating term is derived from (ID.5t : 



+ ~ [Stf^l + 31 (S 2 &1 + ClTo 2 )] + ^ \m 3 D 3 + ^ (Bl&l + Cffi) 
where the S l Q , C l Q , D' Q , £,[ and Tq with I = {1, 2, 3} are given by: 



(E.25) 



^l 1 = & 
gi _ m 

°0 - « 



2^ 



-A, a 



2^ 



So = [7^o] 



T 1 
' 



^?[ 2 = 2i 
« 2 = 2iL 
^2 _ % 



6^ 



-Aa % 



7)2 



6 X 4 



T 2 

-* n 



and 



<??3 _ i?v 



r>3 _ w 



^0 



12f 

^4 



-A 3 1 a 
124 



So = [7^4] 



7 



(?) 



(E.26) 



E.3.2. Baroclinic relation 
We apply (|64|i to I = 2: 



<pA 2 - 6 X V 2 = - 
5 



s 2 + ^£ + Il( r ^ 

rp r dr \ p 



(E.27) 



where T> 2 



rd r Ci 



[rd r + |Q^2 and where Xf^ 2 and J/^ ; 2 have been given in IE.21IE.2*2l . 



E.3.3. Horizontal fluctuation of the molecular weight 

The evolution equation for the mean molecular weight (Paper I, eq. 39) is unchanged by the introduction of the magnetic field: 



dA2 d In fi 
dt dt 



(E.28) 
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E.4. System fori = 4 
E.4.1. Horizontal shear 

We have: 

p- (r 2 Cl 2 ) ~ 2pUr —d r (pr 2 U 2 ) - aU 2 
where the expression for T 2 has been given in jE.9> . 



= -iopv A n 2 + r 2 



E.4.2. Meridional circulation 

In the same way as for I — 2, we recast ( 1771781 in two first order equations as 



U 4 = 



Mg 



1 



' ad 



where 



r 4 = 2 



pC P T) 

7c + 7x ( ? + 



g4 /c,4 + /£,4 fc+ft , xxv , * \ Pm 
— + , , n - (S^A + ^4) + — 

g 4-nGp 4-nGp p 



r „ 20// r / 

-3^U-— — 1 
3 3r \ 



{^14 + (f e er - f e 6 + 6m + (/^ + / eV - <p)A 4 } + y 



and 



yt 4 = i/rS^ - (1 - ^ +Xrm ~ (v +^)A 4 - 

The coefficients related respectively to the centrifugal force and to the Lorentz force are given by: 



fcA = ^d r (r 2 a 4 ) + 20y and fa = ^d r (r 2 ^ 1 ) + 20 



rpo 



where: 



«4 



•|rQ0 2 



fe 4 = £rCX2 2 



and 

^7i;4 



5F {l [( C o £ o ~~ ^o^o) ~~ ( C o £ o ~ ®o?o )] + T ( C o £ o ~ S o^?) + T5T ( C o £ o ~ ^o^) 
We get the relative fluctuation of the effective gravity from i70\ : 



84 



lg 1 / J/r x ;4\ 1 / ^;4 

-j — 7 r & 4 + — — + — 04 + -^- 

^ \ Po / go \ Po 



Po / go 
where 4 obeys the Poisson equation 



d 

H 

dr 



ld 2 r > 20- 

7d72 W-^4- 



47rG dpo— 47rG 
go dr go 



d d / x 

pofl4 + {rp{)b A ) + Xr £ -A + -p ( r ^^;4 j 



(cf. eq.EB- 



We proceed as for / = 2 to calculate the Ohmic heating: 



11 



9^Dg + 



28898 



105 



■(Sg6g+C°f£) 



where the S' , C[ y D l , S' Q and with / = {1,2,3} have been given in the previous section in ( IE.26> . 
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E.4.3. Baroclinic relation 

We apply l|64} to Z = 4. We get: 



<pA 4 - 6^4 = = 



I d I if' 



fl4 + ^ + -T^ 



rp 



r dr 



7 r z -A 



P 



where D4 = ^ [rd r (Q^) ~ 2QQ2] and where Xj^.4 and ifj^A have been given in dE.34t , 
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(E.38) 



E.4.4. Horizontal fluctuation of the molecular weight 

This equation is not affected by the introduction of the magnetic field: 



dA4 d In ytz 
df dr 



Ua 

A 4 = — - —D h A4- 



20, 



(E.39) 



These equations are ready to be implemented in a stellar evolution code, together with the boundary conditions discussed in §7. 



